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ABSTRACT 


Using  the  convenient  second-order  interval  properties  of 
a  two-state  semi-Markov  model  for  a  univariate  point  proc- 
ess, an  automated  technique  for  the  estimation  of  the  param- 
eters in  the  model  was  researched  and  discussed.   The  power 
spectral  density  of  intervals  was  estimated  by  the  period- 
ogram  and  a  Kolmogorov-Smirnov  test  of  fit  was  conducted. 
The  asymtotic  exponential  distribution  and  independence  of 
the  periodogram  points  were  used  to  calculate  an  approximate 
likelihood  function.   A  system  of  equations  was  then  formed 
to  find  the  maximum  likelihood  estimates  of  the  parameters. 
Since  closed-form  solutions  for  the  estimates  could  not  be 
found,  an  iterative  method  to  stabilize  initial  guesses  of 
the  parameter  values  was  attempted  with  only  limited  success. 
Results  on  using  Kolmogorov-Smirnov  type  statistics  and  the 
spectrum  of  intervals  to  test  the  fit  of  stochastic  process 
models  to  data  have  also  been  obtained. 


iua 


TABLE  OF  CONTENTS 

I.  INTRODUCTION  5 

II.  BACKGROUND  ANALYSIS  8 

A.  TWO-STATE  SEMI-MARKOV  MODEL  FOR  UNIVARIATE 
POINT  PROCESS 8 

B.  PARAMETER  ESTIMATION  TECHNIQUE  10 

III.  EXPERIMENTAL  APPROACH  15 

A.  UNIVARIATE  TWO-STATE  SEMI-MARKOV  SIMULATION  —  15 

B.  CALCULATION  OF  THE  PERIODOGRAM 17 

C.  SOLVING  SIMULTANEOUS  EQUATIONS  FOR  THE 

■   PARAMETERS 20 

IV.  RESULTS  AND  CONCLUSIONS 24 

COMPUTER  SUBPROGRAMS  29 

LIST  OF  REFERENCES 38 

INITIAL  DISTRIBUTION  LIST  39 

FORM  DD  1473 40 


ACKNOWLEDGEMENT 

The  author  gratefully  acknowledges  the  cooperation  and 
assistance  of  the  IBM  San  Jose  Research  Laboratory  which 
generously  provided  computer  research  facilities  during  a 
curriculum  sponsored,  six  week,  experience  tour.   In  partic- 
ular, the  support  and  guidance  of  Mr.  G.  S.  Shedler,  IBM 
San  Jose,  and  Professor  P.  A.  W.  Lewis,  Naval  Postgraduate 
School,  were  greatly  appreciated. 


I.   INTRODUCTION 

It  is  in  the  nature  of  the  Operations  Research  approach 
to  the  study  of  problems  to  attempt  the  construction  of  a 
mathematical  model  for  the  problem.   Subclasses  of  mathemat- 
ical models  include  stochastic,  i.e.  utilizing  random  vari- 
ables, and  deterministic  models.   If  a  stochastic  model 
seems  appropriate  and  a  general  model  is  proposed,  it  re- 
mains necessary  to  estimate  parameters  of  the  model  from 
data.   Parameter  estimates,  as  well  as  the  general  form  of 
the  model,  usually  come  from  detailed  analysis  of  data  ob- 
served from  the  problem  or  process  under  investigation. 

Several  techniques  utilizing  observed  data  exist  for  the 
estimation  of  parameters  for  stochastic  models.   Typically 
the  methods  of  moments  or  maximum  likelihood  are  used  and 
usually  yield  estimates  with  some  desirable  properties. 
Methods  such  as  these  frequently  require  the  simultaneous 
solutions  to  a  system  of  equations  in  order  to  find  esti- 
mates.  A  number  of  computer  approximation  routines  have 
been  developed  for  the  solution  of  such  systems,  but  their 
usefulness  seems  limited. 

One  proposed  stochastic  model  provided  the  impetus  for 
this  research.   Lewis  and  Shedler  [1973],  while  studying 
page  reference  patterns  in  a  demand  paged  computer  system, 
formulated  a  univariate  two-state  semi-Markov  model  for  the 
process  of  page  exceptions.   Page  exceptions  occur  because  a 


computer  program  which  is  in  execution  has  been  stored  in 
blocks  of  storage  called  pages.   Some  of  these  pages  must  be 
in  core  storage  for  the  program  to  be  executing,  while  the 
remaining  pages  may  be  located  on  peripheral  storage  de- 
vices.  Following  the  execution  of  each  instruction  a  page 
is  referenced  which  contains  the  next  instruction.   If  this 
referenced  page  is  in  core  storage  execution  continues;  how- 
ever, if  the  referenced  page  is  not  in  core  storage  then 
execution  is  interrupted  and  the  referenced  page  must  be 
read  into  core  storage.   This  type  of  interruption  is  refer- 
red to  as  a  page  exception.   Data  for  this  process  was  gen- 
erated by  counting  the  number  of  page  references  occurring 
between  page  exceptions.   Lewis  and  Shedler  [1973]  discussed 
their  procedure  for  estimating  parameters  which  they  de- 
scribed as  an  ad  hoc  method,  and  concluded  that  there  was  a 
need  to  formalize  the  parameter  estimation  procedure. 

The  purpose  of  the  research  in  this  thesis  was  to  uti- 
lize the  convenient  second-order  interval  properties  of  a 
univariate  two-state  semi-Markov  process  to  produce  an  auto- 
mated, computer  programed,  technique  for  the  estimation  of 
parameters  for  the  model.   This  was  desirable  because  the  ad 
hoc  method  used  by  Lewis  and  Shedler  [1973]  was  very  time- 
consuming  and  there  exists  a  considerable  body  of  page  ex- 
ception data  which  it  is  desired  to  analyze.   The  basic 
procedure  was  to  calculate  an  estimate  for  the  power  spec- 
tral density  of  the  process,  namely  the  periodogram,  and 


utilize  an  approximate  method  of  maximum  likelihood  to  esti- 
mate the  parameters. 

It  will  be  seen  that  the  proposed  procedure  did  not  work 
as  well  as  hoped,  but  the  problems  which  arose  pointed  up 
other  possible  attacks  on  the  problem.   It  should  also  be 
noted  that  model  fitting  and  parameter  estimation  for  these 
point  processes  is  almost  a  completely  open  field. 


II.   BACKGROUND  ANALYSIS 

A.   TWO-STATE  SEMI-MARKOV  MODEL  FOR  UNIVARIATE  POINT  PROCESS 

Excellent  discussions  of  this  model  can  be  found  in  Cox 
and  Lewis  [1966, Ch. 7]  and  Lewis  and  Shedler  [1973],   Those 
discussions  are  summarized  here  for  continuity  of  exposition. 

Let  the  sequence  of  random  variables  {X- ,i=l, . . . ,N}  be 
interevent  times,  i.e.  X-  is  the  interevent  time  between 
event  (i-1)  and  event  (i) .   In  order  that  a  discussion  of 
equilibrium  distributions  may  be  avoided  it  was  assumed  that 
a  hypothetical  event  has  occurred  at  time  zero,  so  that  Xj, 
the  interval  between  time  zero  and  the  first  event,  is  an 
observation  from  the  same  process  as  the  remainder  of  the 
sequence,  i.e.  there  is  no  length-biased  sampling  [Cox  and 
Lewis  1966, Ch. 4]  included. 

Now  suppose  there  are  two  types  of  intervals  but  that 
the  interval  type  is  not  observable,  i.e.  a  univariate  point 
process.   The  two  interval  types  have  probability  mass  func- 
tions (p.m.f.)  Pi (x)  and  p^ (x) ,  respectively,  with  transi- 
tions between  types  described  by  a  two-state  Markov  chain 
with  matrix 

A  =   /  °'-     l-« 


That  is,  given  that  X.  has  p.m.f.  Pj (x)  then  X.    has  p.m.f. 
P2 (x)  with  probability  Oj  and  p.m.f.  PjCx)  with  probability 
l-ttj,  independent  of  the  history  of  previous  intervals,  etc. 


The  vector  of  steady-state  probabilities  tt  =  (tTj  tTj) 
associated  with  the  transition  matrix  A  results  from  the 
solution  of  the  matrix  equation  tt  =  jtA  and  it  follows  that 
1-a  1-a 

2-ai-a2  2-ai-a2 

If  \i^,a^^   and  1^2' "^2^  ^^®  ^^^   mean  and  variance  for  intervals 
with  p.m.f.  Pi(x)and  P2 (x) ,  respectively,  the  steady-state 
marginal  results  for  intervals  between  events  in  the  univar- 
iate process,  i.e.  interval  type  not  known,  are  as  follows: 
P(x)  =  TTiPi(x)  +  7r2P2(x)  , 

y     =    E(X)      =     TTlUl     +     ■IT2y2      , 

a^   =  var(X)  =  t^iOi^   +   -nzOz^    +   tti1T2  (Ui-ya)  ^  • 

The  serial  correlation  coefficients  of  lag  k,  p,  ,  for  the 

k   9 
intervals  are  of  the  form  mB  /a^  where,  for  k=l,2,..., 

m  =  TTiirz  (yi-ya)  ^        6  =  Uj  +  ttj  -1  • 
From  these  coefficients  the  positive  portion  of  the  power 
spectral  density  may  be  computed, 

P  {(jj  )  =  TT   (1  +  2   Z   p,  cos  kui)  . 
+  n  k=l   ^      "^ 

The  closed-form  solution  to  the  infinite  series  is  given  by 

Jolley  [1961,  series  #545]  yielding 

o^  mB    (cos  03  )  -  B 

P+{u)j^)  =  —  [1  +  2  —  { }].    (1) 

IT         o^      l+B^-2Bcos  00^ 

The  beneficial  feature  of  the  power  spectral  density  for 
this  model  is  that  it  only  depends  upon  the  mean  and  vari- 
ance of  each  of  the  two  probability  distributions  and  not 
on  the  complete  distributions,  and  is  thus  fairly  robust. 


The  count  spectrum  [Cox  and  Lewis  1966,  Ch.4]  on  the  other 
hand,  depends  on  the  complete  distributions. 

B.   PARAMETER  ESTIMATION  TECHNIQUE 

Lewis  and  Shedler  [1973]  used  a  modified  method  of  mo- 
ments approach  in  order  to  estimate  the  parameters  in  their 
model.   The  standard  method  of  moments  procedure  for  param- 
eter estimation  is  to  calculate  theoretical  moments  in  terms 
of  the  unknown  parameters  and  equate  them  to  empirical  ob- 
servations of  these  moments.   An  alternative  to  this  method 
is  the  method  of  maximum  likelihood.   In  this  method  param- 
eter values  are  selected  which  maximize  the  joint  probabil- 
ity density  of  the  observed  data.   To  accomplish  this  there 
is  a  need  for  some  distributional  assumptions.   However, 
even  for  a  simple  model  such  as  the  univariate  two-state 
semi-Markov  model  discussed  here  it  is  not  possible  to  write 
down  the  joint  density  of  the  intervals.   Thus  the  following 
approximate  technique  was  proposed  and  tried. 

It  is  known.  Cox  and  Lewis  [1966, Ch. 5],  that  an  estimate 

of  the  power  spectral  density  P  (oj  )  at  oj  ,  the  periodogram 

I(n),  is  in  general  asymtotically  exponentially  distributed 

[Olshen  1967] .   The  periodogram  is  an  unbiased  estimate, 

i.e.  E[I(n)]=P((jj  );  however,  it  is  not  a  consistent  estimate 

n 

since  the  variance  of  the  exponential  distribution  is  equal 
to  the  square  of  the  mean,  i.e.  the  variance  does  not  de- 
crease with  increased  sample  sizes.   Moreover,  for  nj  not 
equal  to  nj  the  periodogram  points  I(ni)  and  iCnj)  are 
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asymtotically  independent.   Thus  for  finite  sample  size  N  an 
approximate  likelihood  function  may  be  written  by  assuming 
the  periodogram  points  are  independent  with  exponential  dis- 
tributions having  mean  value  P(cjo  ).   This  is  the  technique 
explored  in  this  thesis. 

The  definition  and  development  of  the  periodogram  re- 
quires the  finite- Fourier  transform. 

The  finite  Fourier  transform  was  discussed  by  Cooley, 
Lewis  and  Welch  [to  be  published  in  1974].   Let  {Y(j), 
j=0,...,N-l}  be  a  sequence  of  N  real  numbers.   The  finite 
Fourier  transform  of  Y(j)  is  then 

N-1      _2TTini 
a(n)  =  i  Z  Y(j)e   ^   ^  n=0,...,N-l. 
N  j=0 

This  sequence  of  complex  numbers  may  also  be  written  in  the 

form 

.   N-1  N-1 

a(n)  =  -   2  Y(j)cos(jw  )-i  Z  Y(j)sin(j{jo  )  , 
N   j=0  ""        j=0  ^ 

where  o)  =  2vn,    n=0,l,  .  .  .  ,N-1.   The  periodogram  I  (n)  is  then 

~N 

N|a(n)  1^ 
I  (n)  =   2Tr     ,  n=0,  .  . .  , 

or 

N-1  N-1 

I(n)  =  {  Z  Y(j)cos(ja3„)}^  +  {  E  Y(j)sin{ja)  )}^ 

JizO 1=0 1 

27TN 

The  periodogram  is  an  even,  periodic,  function  and  hence  has 

only  [N/2]+l  distinct  values,  where  [N/2]  is  the  integer 

part  of  N/2.   Hereafter  in  the  discussion  N  will  refer  to  an 

even  integer.   Let  I,  (n)=2l(n)  be  the  estimate  for  P,  (w  ). 

+  +   n 
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It  is  easily  seen  that  I  (0)  is  proportional  to  the 
square  of  the  arithmetic  average  of  the  observed  data;  thus 
no  new  information  is  obtained  from  I  (0) .   Since  N/2  is  an 
integer  ^f^z-y^'^   and  I  (N/2)  is  proportional  to  the  square  of 
an  alternating  summation  of  the  data.   Both  I  (0)  and 
I  (N/2)  were  ignored  in  what  follows,  thus  leaving  (N/2)-l 
periodogram  points.   It  should  be  added  that  these  two  peri- 
odogram  points,  suitably  normalized,  have  asymtotic  x^  dis- 
tributions with  one  degree  of  freedom  and  not  an  exponential 
distribution. 

Now  there  is  sufficient  information  to  begin  the  approx- 
imate maximum  likelihood  search  for  parameter  estimates. 
The  parameters  of  this  model  that  need  estimation  are  the 
mean  and  variance  of  each  marginal  distribution  and  the  two 
transition  probabilities  aj  and  o-z.      As  a  vector  these  pa- 
rameters will  be  labeled  9^=  (y  j  ,ai  ^  ,y2 'Cfa  ^ '"^i 'C'a)  ^'^'^   indi- 
vidually, to  simplify  notation,  as  6  ■ , j=l , 2 , 3 , 4 , 5 , 6 ,  to 
stand  for  the  parameter  as  an  element  of  the  vector  £. 

The  approximate  likelihood  function  can  be  written  as 

(N/2)-l     1       -I+(n)/P,  (aj^;e) 
L(£)  =    n      P+(aj  ;0_)   e  , 

n=l 

which  is  equivalent  to 

(N/2)-l  I  (^) 
/(N/2)-l     1     \  -I,  V^  ^ 

L(e)  =  (  n    pTT^vIT)  ^     P+f'^n'-i)  • 

\  n=l  / 

A  more  simple  function  to  work  with,  which  has  the  same  max- 
imum as  L  ( e_)  ,  is"  the  log  likelihood  function  LL(e_)=ln  L  (e_)  , 
where  In  symbolizes  the  natural  logarithm. 
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Then 

(N/2)-l  (N/2)-l  ^  ,  ^ 

LL(e_)  =  -  E  ln{P  (oj  ;6)}   -    I  i-+(n) 

n=l        "^  "      n=l    P+(%;i)  . 

In  the  typical  mathematical  approach  to  finding  an  un- 
constrained maximum  of  a  function,  it  is  a  necessary  condi- 
tion that  all  of  the  first-order  partial  derivatives  of  the 
function,  with  respect  to  the  unknown  parameters,  be  equal 
to  zero,  i.e.  that 

0  =  9LL  ( e )  =  LL^  =  E      I+(n)-P4.(a)ri;B)  p.,j=i,...,6    (2) 
aej       -•   n=l     P+(a)j^;i)       ^ 

where  P-  and  LL .  are  the  first-order  partial  derivatives  of 
P  (u3  ;£)  and  LL  (£)  ,  respectively,  with  respect  to  parameter 
6..   This  process  results  in  six  equations  and  six  unknown 
parameters.   Parameter  estimates  are  found  by  simultaneously 
solving  the  system  of  equations  for  each  of  the  parameters, 
although  this  may  not  yield  a  unique  maximum.   If  the  system 
is  of  a  simple  form  it  may  be  possible  to  get  at  least  a  few 
closed-form  solutions  which  will  reduce  the  size  of  the  sys- 
tem. 

Once  the  parameter  estimates  have  been  found  it  is  nec- 
essary to  show  that  a  maximum  has  been  achieved.   A  suffi- 
cient condition  for  a  maximum  is  that  the  matrix  of  second- 
order  partial  derivatives  be  negative  definite.   The  final 
phase  in  this  approximate  likelihood  estimation  process  is 
to  verify  the  predictability  of  the  model.   The  verification 
may  be  done,  using  the  estimated  parameters,  by  calculating 
other  theoretical  properties  of  the  model,  such  as  the  spec- 
trum of  counts  discussed  by  Cox  and  Lewis  [1966,  Ch.4], 
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which  may  then  be  compared  with  the  corresponding  empirical 
properties  of  the  data.   Note  that  the  utility  of  the  spec- 
trum of  intervals  in  the  approximate  likelihood  analysis  is 
that  it  does  not  depend  on  the  complete  distributional  form 
for  pi (x)  and  pa (x)  while  the  spectrum  of  counts  does.   It 
will  be  seen  later,  however,  that  this  independence  leads  to 
ill-conditioning  in  the  solution  of  the  maximum  likelihood 
equations. 
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III.   EXPERIMENTAL  APPROACH 

The  original  data,  analyzed  by  Lewis  and  Shedler  [1973] , 
was  not  available  for  this  research.   In  view  of  this  fact 
and  since  the  purpose  of  the  research  was  to  evaluate  the 
effectiveness  of  the  previously  described  technique  for  pa- 
rameter estimation  it  was  felt  that  data  observed  from  a 
model  with  known  parameters  would  better  aid  the  evaluation 
process.   With  this  in  mind  a  simulation  of  the  model  de- 
scribed by  Lewis  and  Shedler  [1973]  was  constructed  for  the 
purpose  of  generating  such  data. 

A.   UNIVARIATE  TWO-STATE  SEMI-MARKOV  SIMULATION 

The  simulation,  as  well  as  the  model,  was  subdivided 
into  three  major  subsections.   The  state  transition  matrix  A 
was  one  subsection  and  the  two  distributions  for  intervals 
were  the  remaining  two  subsections.   Lewis  and  Shedler  [1973] 
postulated  a  geometric  distribution  for  the  long  intervals 
and  a  negative  binomial  distribution  for  the  shorter  inter- 
vals.  The  parameters  used  for  the  simulation  were  those 
calculated  by  Lewis  and  Shedler  [1973] . 

A  Monte  Carlo  simulation,  such  as  this,  required  a 
pseudo-random  number  generator  with  favorable  serial  corre- 
lation properties.   Learmonth  and  Lewis  [1973]  discussed 
such  a  generator  called  SRAND.   SRAND  returns  an  observation 
from  a  standardized  uniform  distribution  on  the  interval 
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(0,1).   SRAND  is  a  multiplicative  generator  with  a  multi- 
plier of  (7=)  and  a  modulus  of  (23^-1). 

The  geometric  distribution  is  of  the  form 

Pi(x)  =  pj^"^  (1-Pi)/  0<pj<l;  x=l,2,..., 
with  a  mean  yj=l/(l-pj)  and  variance  Cj ^=Pj/ (1-pj ) ^ .   Uti- 
lizing the  survivor  function  of  the  geometric  distribution, 
i.e.  prob{X>x}=p j^,  x=l,2,...  ,  a  generator  of  geometric 
variates  was  obtained.   It  was  of  the  form 

X  =  r{ln(R)/ln(Pi)}  , 
where  R  was  an  observation  from  SRA.ND  and  the  symbol  [  {h} 
signified  the  smallest  integer  greater  than  or  equal  to  b. 

The  negative  binomial  distribution  is  of  the  form 

PaCx)  =  f""  "]   P,^~^  (l-P,)   / 


=(r;) 


0<P2<1/  k>0,  x=l,2,...  ,  with  mean  y2=l+{kp2/ (l-Pa ) }  and 
variance  a2^=kp2/ (I-P2 ) ^ •   Let  X|X,  denoting  X  given  a  fixed 
value  of  X,  be  distributed  as  a  Poisson  random  variable  with 
parameter  X.   Now  let  X  have  a  gamma  distribution  with  pa- 
rameters k  and  r], 

k  k-i  -nx 

f  (X)  =  -2 ^ ,  k>0;  X>0;  n>0. 

r(k) 

It  can  be  shov/n  using  generating  functions  that  the  uncondi- 
tioned X  has  a  negative  binomial  distribution  with  parame- 
ters k  and  p2=l/(l+n). 

To  calculate  a  gamma  variate  with  a  parameter  k,  a  posi- 
tive real  number,  it  was  necessary  to  employ  Johnk ' s  tech- 
nique [1964]  for  generating  variates  with  the  fractional 
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part  of  k.  Let  k  be  the  integer  part  of  k,  if  k>^l,  or  zero 
if  k<l,  and  let  R"  be  the  fractional  part  of  k.  The  sum,  Xj, 
of  k  exponentially  distributed  random  variables  with  param- 
eter n  has  a  gamma  distribution  with  parameters  k  and  n-  In 
Johnk's  technique  let  Uj,  Uj  and  U3  be  independent  and  iden- 
tically distributed  observations  from  a  uniform  distribution 
on  interval  (0,1). such  that 

Y  =  u//^+  U,^/(^-'^^  <1  . 
If  Y>^1  new  observations  for  Uj  and  U2  should  be  obtained. 
Then  for  Z=Ui^/^/Y  and  E=-ln  U3,  X2=(2xE)/n  has  a  gamma  dis- 
tribution with  parameters  k  and  n.   Finally  A.=Xj+X2  has  the 
required  gamma  distribution  with  parameters  k  and  n . 

The  generation  of  Poisson  random  variates  with  parameter 
X  was  accomplished  by  letting  X  be  equal  to  the  largest  in- 
teger n  such  that,  for  a  sequence  of  independent  identically 
distributed  uniform  random  variates  (U. )  from  the  interval 


(0,1), 


U  xu  X. . .xu  >e"^  . 
12  n 


If  Ui_<e~   then  X=0.   X  is  then  distributed  as  a  Poisson  var- 
iate  with  parameter  X. 


B.   CALCULATION  OF  THE  PERIODOGRAM 

The  finite  Fourier  transform  discussed  in  section  II. B 
above  requires  on  the  order  of  N^  complex  operation  pairs, 
i.e.  a  multiplication  and  an  addition.   For  large  N  this 
can  be  very  costly  in  terms  of  calculation  time.   Cooley, 
Lewis  and  Welch  [1970]  discussed  the  use  of  a  fast  Fourier 
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transform  algorithm  which  only  requires  on  the  order  of 

N(ri+...+r  )  complex  operation  pairs  where  N= (r ^x. . . xr  ) , 

i.e.  r   is  a  factor  of  N.   The  International  Mathematical 
m 

and  Statistical  Library  [1973  revision]  contains  a  computer 
subroutine,  FFTR,  which  computes  the  fast  Fourier  transform 
of  a  real  data  sequence.   For  N=820,  as  in  this  research, 
the  fast  Fourier  transform  algorithm  used  only  six  percent 
of  the  number  of  complex  operation  pairs  required  by  the 
straight- foirward  calculation  method.   Thus  a  significant 
savings  in  computer  operating  time  was  realized. 

Utilizing  previously  described  equations  the  periodogram 
I+(n)  was  computed  and  then  used  in  a  test  of  fit  to  the 
power  spectral  density  P^(cOj^).   Cox  and  Lewis  [1966,  Ch.6] 
described  a  test  based  on  the  uniform  distribution.   While 
the  periodogram  has,  asymtotically,  an  exponential  distri- 
bution with  mean  P,(u  ),  the  quantity  I^(n)/P  (u  )  has  an 
exponential  distribution  with  mean  one.   This  is  true  for 
each  of  the  (N/2)-l  periodogram  points.   If  all  (N/2)-l  of 
these  quantities  are  summed  the  total  gives  an  interval  of 
length  over  which  there  are  (N/2) -2  points  dispersed.   The 
intervals  between  these  points  are  each,  hypothetically,  an 
observation  from  a  unit  exponential  distribution,  i.e.  the 
points  form  a  Poisson  process.   It  is  a  well  known  fact  of 
the  Poisson  process  that  given  M  points  are  in  an  interval 
the  M  points  are  dispersed  uniformly  over  the  interval. 
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Thus,  the  sequence  tU/-*,  i=l, . . . , (N/2)-2},  where 

i 
j^I^  I+(n)/P+(a)n) 


U 


(i)  ~(N/2)-l 


I.         I  (n)/P_^{(j  ) 
n=l    +     +  n 


are  uniform  order  statistics.   The  empirical  cumulative  dis- 
tribution function  for  these  quantities  was  then  compared 
with  the  uniform  cumulative  distribution  function  using  the 
Kolmogorov-Smirnov  test.   The  null  hypothesis  is  that  the 
sequence  {U/j_\}  is  formed  of  uniform  order  statistics,  while 
the  alternative  hypothesis  remains  unspecified.   Lilliefors 
[1969]  found  that  the  critical  values  of  the  K-S  test  are 
too  conservative  when  testing  using  exponential  distribu- 
tions where  the  mean  has  been  estimated,  as  in  this  case. 
Too  conservative  means  that  the  listed  critical  value  for  a 
level  of  significance  a  has  actually  a  level  of  significance 
less  than  a.   If  the  above  test,  with  modified  percentage 
points,  accepts  the  null  hypothesis  then  the  assumption  of  a 
semi -Markov  model  for  the  data  has  been  justified. 

In  order  to  test  the  periodogram  it  was  necessary  to 
know  P^(a)j^).   As  discussed  earlier  the  correlation  coeffi- 
cient of  lag  k,pj^,  is  pj^=m3  /cJ^  for  this  model.   Let  y(0), 
Y(1)  and  y(2)  be  estimates  of  the  variance  and  covariances 
of  lags  one  and  two,  respectively,  for  the  intervals.   Then 

Yd)  =  o^Pj  =  m3 
and 

Y(2)  =  a^P2  =  mB^   . 
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Solving  simultaneously  for  m  and  3/  the  estimates  of  m  and  3 

are 

I   =  y(2)       and       m  =  y' (1) 

ytlT  yUT      . 

From  (1)  an  estimate  for  P.  (%)  ^^s 

1^  ^    (cos  0)  )  -R 

P+(cJn)  =  ^  ^^(0)  +  2mB } 

1+3^-2b  cos  Wn 

These  estimates  were  then  used  in  the  computations  for  the 
sequence  ^^U  . . «  }. 

C.   SOLVING  SIMULTANEOUS  EQUATIONS  FOR  THE  PARAMETERS 

The  system  of  equations  defined  by  (2)  and  (1)  was  ex- 
tremely complex,  with  no  hope  of  finding  a  closed- form  solu- 
tion for  any  of  the  parameters.   The  system  was  reduced, 
however,  by  noting  from  the  geometric  distribution  assump- 
tion that  the  variance  for  the  long  intervals  was  a  function 
only  of  the  parameter  pi,  which  also  was  the  only  parameter 
in  the  mean.   It  was  a  simple  matter  then  to  find  the  vari- 
ance as  a  function  of  the  mean  which  then  reduced  the  system 
to  only  five  unknown  parameters.   The  system  was  still  com- 
plex and  required  some  iterative  method  for  solution. 

Rao  [1965]  suggested  an  iterative  method  which  he  called 
the  Method  of  Scoring.   He  called  LL.  the  jth  efficient 
score.   The  approach  for  this  method  was  to  assume  some  ini- 
tial trial  solution.   Using  a  first-degree  Taylor's  expan- 
sion of  the  efficient  scores  about  the  trial  solution,  a 
system  of  linear  equations  was  derived  from  which  an  addi- 
tive correction  to  the  trial  solution  was  found.   The 
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iterations  were  repeated  until  the  additive  corrections  be- 
came negligible. 

Specifically,  let  Q ^^ , . . .  ,Q ^°   be  the  trial  values  for 
the  unknown  parameters.   From  (2) 

3LL (i)      5  3^LL(a) 

LL.    ae."  +  .s  (9i-ei°)3e.o90  0  ,  j=i,...,5  , 

J         J       1=1  J     -^ 

where  3LL(£)/3e  .  °==S  . "  the  first-order  partial  derivative  of 
LL(£)  with  respect  to  0.,  evaluated  at  Q.°.      Let  66 j= (6 .-6 . ") 
and  also  let  3^LL(e_)/ (39  j  °3ei°)=Tji° .   Then 

5 

-  E  T.i''59i=  S.°,    j=l,...,5 
i=l  -"^       =• 

was  a  system  of  linear  equations  with  five  unknowns.   In 

matrix  notation  the  system  had  the  form  -T66=S.   Finally, 

the  additive  corrections  were  obtained  from  the  equation 

6_9=-T~^S  where  T~^  was  the  inverse  of  the  matrix  T,  assuming 

T  was  nonsingular.   The  new  trial  solution  then  became 

e_^  =  g^+e^. 

Rao  [1965]  explained  that  the  variance  of  the  final  es- 
timate 6^^  of  9.  was  approximated  by  the  jth  diagonal  ele- 
ment of  the  matrix  (-T~M.   Recalling  that  the  matrix  of 
second-order  partial  derivatives  of  LL  (e_)  should  be  negative 
definite,  then  -T  and  (-T~M  were  both  positive  definite. 

In  order  to  apply  the  method  of  scoring  it  was  necessary 
to  determine  initial  estimates  of  the  parameters.   The  mean 
and  variance  of  the  intervals  and  the  parameters  m  and  3  all 
have  been  estimated.   Utilizing  the  marginal  properties  of 
the  model  and  the  method  of  moments  a  system  of  four 
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equations  was  developed  which  was  of  the  form 
X  =  (1-S2)li  +  (l-iillz 


(l-g) 

2 


Y{0)  =  a-^z)   (Pi^-Ui)  +  (l-aQg;'  +  m   ; 


m  =  (1-ai)  (1-52)  (^1-^2)^ 

where  X  was  the  estimate  for  the  mean  of  the  intervals  and 
the  quantity  (yj^-yi)  was  the  estimate  for  a ^^    after  making 
the  geometric  distribution  assumption.   From  this  system 
initial  estimates  for  four  parameters,  as  functions  of  the 
fifth  parameter,  were  found  and  had  the  form 
Xyi  -  P  -  m 


V2 

(Pr 

-X) 

• 

= 

(1- 

e)x 

- 

(y,- 

■Byi) 

"l 

(yr 

-Ui) 

1 

02 

= 

3+1 

-Si 

t 

(1-6)  (Y(O)-m)  -  (l-a^)  (y/-yi) 
(1-ai) 
It  only  remained  then  to  estimate  parameter  y j . 

Lewis  and  Shedler  [1973]  explained  that  their  estimate 
of  Pi  involved  an  eyeball  judgement  of  where  linearity  began 
in  the  tail  of  the  log  survivor  function  of  the  data.  This 
linearity  in  the  tail  led  to  the  postulate  of  the  geometric 
distribution  for  the  long  intervals.  Since  only  an  initial 
estimate  was  needed,  their  method  of  estimating  y i  was  uti- 
lized again.   The  value  of  the  interval  where  linearity 
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began  was  subtracted  from  each  greater  interval  and  the 
arithmetic  average  of  these  intervals  was  then  taken  as  the 
estimate  of  ui.   Now  all  parameters  had  been  initially  esti- 
mated and  the  iterative  method  of  scoring  was  applied  to 
stabilize  these  estimates. 
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IV.   RESULTS  AND  CONCLUSIONS 

As  the  research  for  this  thesis  progressed  two  areas 
developed  results  which  need  to  be  discussed.   The  first  of 
these  was  the  test  for  justification  of  the  exponential  dis- 
tribution assumption  for  the  periodogram  points.   A  subrou- 
tine called  KSTEST  was  written  to  conduct  this  test.   Part 
of  the  output  of  this  subroutine  was  the  Kolmogorov-Smirnov 
statistic  which  then  was  compared  to  a  critical  value  from 
the  distribution  proposed  by  Lilliefors  [1969]  for  the  case 
where  the  mean  of  the  exponential  distribution  had  to  be  es- 
timated.  The  0.99  quantile  of  that  distribution,  i.e.  a  one 
percent  level  of  significance,  was  1.25.   It  was  noted  that, 
at  this  level,  of  four  thousand  trials  made  approximately 
six  percent  were  rejected  as  not  having  produced  periodo- 
grams  from  a  semi-Markov  model. 

In  addition  to  testing  the  hypothesis  for  each  simula- 
tion another  benefit  was  received.   Since  the  testing  did 
not  strictly  conform  to  that  discussed  by  Lilliefors  [1969] , 
because  each  periodogram  point  had  a  different  mean,  it  was 
felt  that,  for  this  case,  quantiles  of  the  distribution 
should  be  estimated.   The  four  thousand  data  points  of  the 
statistic  were  obtained  from  four  computer  runs  each  contain- 
ing one  thousand  simulations.   For  each  run,  the  data  was 
divided  into  ten  sections  in  serial  order,  i.e.  the  first 
one  hundred  points  were  the  first  section,  etc.   The  ele- 
ments of  each  section  were  ordered  and  the  0.80,  0.85,  0.90, 
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0.95  and  0.99  empirical  quantiles  were  observed.   This  re- 
sulted in  ten  observations  for  each  quantile  from  which  a 
mean  and  variance  were  estimated.   Lastly,  the  entire  data 
for  the  run  was  ordered  and  the  five  quantiles  were  observed, 
Thus,  for  each  run,  each  of  the  five  empirical  quantiles  had 
been  observed  and  had  an  estimated  mean  and  variance.   Fin- 
ally a  mean  of  the  four  overall  observations  for  each  quan- 
tile and  a  mean  of  the  section  means  were  computed.   The 
results  are  shown  in  Table  I. 

Lilliefors  [1967]  discussed  the  Kolmogorov-Smirnov  test 
for  normal  data  and  calculated  numerically  the  quantiles  of 
the  test  statistic  for  the  case  where  the  mean  and  variance 
of  the  normal  distribution  must  be  estimated.   These  quan- 
tiles are  included  for  comparison. 

The  second  of  the  two  significant  areas  was  the  estima- 
tion of  parameters.   The  subroutine  ESTIM8  was  written,  in 
double  precision,  to  utilize  the  method  of  scoring  for  pa- 
rameter estimation.   As  a  result  of  the  use  of  the  subrou- 
tine several  potential  hazards  to  the  proposed  technique 
became  visible. 

The  first  of  these  hazards  was  the  disparity  between 
magnitudes  of  the  five  unknown  parameters.   Three  of  these 
are  means  and  variances  and  the  other  two  are  probabilities, 
which  are  always  less  than  or  equal  to  one.   This  problem 
became  apparent  when  the  magnitudes  of  the  scores  and  the 
elements  in  the  matrix  of  second-order  partial  derivatives 
were  seen.   An  attempt  to  correct  this  problem  was  made  by 
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TABLE  I 


Source 


0.20 


Level  of  Significance 
0.15    0.10    0.05 


0.01 


Usual  quantiles 

1.07 

1.14 

1.22 

1.36 

1.63 

Lilliefors  quantiles 

0.86 

0.91 

0.96 

1.06 

1.25 

Run  1 

0.74 

0,79 

0.96 

1.66 

8.56 

Mean 

0.73 

0.80 

0.92 

1.59 

6.04 

Variance 

0.002 

0.005 

0.01 

0.31 

13.24 

Run  2 

0.75 

0.85 

1.03 

2.05 

6.66 

Mean 

0.75 

0.85 

1.03 

1.86 

7.11 

Variance 

0.003 

0.01 

0.04 

0.45 

25.35 

Run  3 

0.74 

0.80 

0.90 

1.31 

6.79 

Mean 

0.73 

0.79 

0.91 

1.80 

4.90 

Variance 

0.001 

0.002 

0.01 

1.31 

9.87 

Run  4 

0.73 

0.82 

0.96 

1.72 

8.45 

Mean 

0.74 

0.83 

0.95 

1.64 

5.90 

Variance 

0.01 

0.02 

0.03 

0.62 

12.51 

Mean  of  Runs 

0.74 

0.82 

0.96 

1.69 

7.62 

Variance  of  Runs 

0.00 

0.00 

0.001 

0.03 

0.27 

Mean  of  Means 

0.74 

0.82 

0.95 

1.72 

5.99 

Variance  of  Means 

0.00 

0.00 

0.001 

0.01 

0.21 

Lilliefors  normal 

quantiles 

0.736 

0.768 

0.805 

0.886 

1.031 
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dividing  the  three  large  parameters  by  the  overall  mean  of 
the  intervals.   This  also  favorably  affected  the  partial 
derivatives  involving  these  parameters.   The  desired  effect 
was  achieved  in  that  the  gap  between  magnitudes  was  nar- 
rowed; however,  the  parameter  estimates  that  resulted  from 
this  modification  were  only  about  one  percent  different  from 
the  parameters  achieved  earlier,  so  apparently  the  disparity 
created  no  significant  problem. 

The  second  of  these  problems  was  that  the  final  param- 
eter estimates,  overall,  appeared  to  have  little  relation- 
ship to  the  marginal  parameter  values  from  which  the  data 
was  generated.   Similarly,  parameter  estimates  for  two  sets 
of  data  differed  greatly  in  magnitude  and  at  times  in  sign, 
even  when  the  periodogram  was  accepted  as  a  close  fit  to  the 
power  spectral  density.   Differences  in  sign  were  extremely 
disturbing  since  all  of  the  parameters  were  expected  to  be 
greater  than  zero. 

A  third  problem,  related  to  the  second,  was  that  the  re- 
sults failed,  numerically,  to  establish  that  the  matrix  of 
second-order  partial  derivatives  was  negative  definite. 
Similarly  the  negative  inverse  of  that  matrix  could  not  be 
shown  to  be  positive  definite.   This  problem  indicated  that 
either  a  maximum  had  not  been  achieved,  even  though  a  cut- 
off criterion  of  10"'"  was  used  to  test  for  convergence,  or 
that  due  to  round-off  error  the  properties  of  a  maximum 
could  not  be  detected.   With  a  smaller  cut-off  criterion  the 
process  would  not  converge  and  had  to  be  terminated. 
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Finally,  in  a  few  instances  when  four  of  the  five  un- 
known parameters  appeared  close  to  the  simulation  parameters 
the  initial  value  of  the  fifth  parameter  was  changed  and  the 
subroutine  was  restarted.   The  parameters  would  again  con- 
verge; however,  the  final  values  in  some  cases  changed  dras- 
tically, even  to  the  point  of  changing  sign. 

Some  of  these • problems  may  have  been  caused  by  an  ill- 
conditioned  system  of  equations,  while  others  might  be  due 
to  the  lack  of  a  powerful  iterative  technique  for  the  solu- 
tion of  a  system  of  equations  that  has,  perhaps,  poor  initial 
estimates.   In  any  case  it  should  be  clear  that  the  use  of 
second-order  properties  of  a  model  might  simplify  or  at 
least  aid  the  parameter  estimation  process.   One  proposed 
modification  to  the  technique  discussed  in  this  thesis  was 
to  use  a  mixture  of  the  method  of  moments  approach  on  the 
marginal  distribution  of  the  intervals  and  the  maximum  like- 
lihood approach  on  the  second-order  properties  to  estimate 
parameters. 

In  conclusion  it  should  be  recalled  that  model  fitting 
and  parameter  estimation  for  univariate  point  processes  is 
almost  a  completely  open  field  and  that  attempts,  even  un- 
successful ones,  are  needed  in  order  to  break-through  the 
barrier  of  inadequate  methodology. 
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COMPUTER  SUBPROGRAMS 


C                       SUBROUTINE  MODEL  C 

c  c 

C    A.   PURPOSE:  C 

C        THIS  SUBROUTINE  CONTROLS  THE  SIMULATION  OF  A  C 
C        UNIVARIATE  TWO-STATE  SEMI-MARKCV  POINT  PROCESS.     C 

C        THE  STATE  ONE  INTERVALS  HAVE  A  GECMETRIC  C 

C        DISTRIBUTION  AND  THE  STATE  TWO  INTERVALS  HAVE  C 

C        A  NEGATIVE  BINOMIAL  DISTRIBUTION.  C 

C  C 

C    e.   USAGE:  C 

C        1.   ARGUMENTS:  C 

C  C 
C            X  -  OUTPUT  VECTOR  OF  INTEREVENT  TIMES  (REAL*8)  C 

C  C 

C            SIZE  -  INPUT  LENGTH  OF  VECTOR  X  (INTEGER)  C 

C  C 

C            IX  -  INPUT  RANDOM  NUMBER  SEED  (INTEGER)  C 

C  C 

C             DISTIM  -  INPUT  MEAN  OF  THE  GECMETRIC  C 

C                      DISTRIBUTION  (REAL*8)  C 

c  c 

C  DIST2M  -  INPUT  MEAN  OF  THE  NEGATIVE  BINOMIAL    C 

C                      DISTRIBUTION  (REAL*8)  C 

c  c 

C            DIST2K  -  INPUT  PARAMETER  K  IN  THE  NEGATIVE  C 

C                      BINOMIAL  DISTRIBUTION  (REAL*8)  C 

C  C 

C            Al  -  INPUT  PARAMETER  IN  THE  TRANSITION  C 

C                      MATRIX  (REAL*8)  C 

C  C 

C            A2  -  INPUT  PARAMETER  IN  THE  TRANSITION  C 

C                      MATRIX  (REAL*8)  C 

C  C 

C        2.   REQUIRED  SUBPROGRAMS:  C 

C            INTEGER  FUNCTION  GEOMET  C 

C            INTEGER  FUNCTION  NEGBIN  C 

C            SUBROUTINE  OVFLOW  (NPS  ROUTINE)  C 

C            SUBROUTINE  SRAND   (NPS  ROUTINE)  C 

C            SUBROUTINE  SNORM   (NPS  ROUTINE)  C 

C            SUBROUTINE  SEXPON  (NPS  ROUTINE)  C 

c    , c 

c  c 

SUBROUTINE  MODEL ( X , S I ZE t I  X, DI ST IM ,DI ST2M , D I ST2K, Al , A2) 

IMPLICIT  REAL*8  (A-H,K,0-Z) 

REAL*4  R 

INTEGER*4  GEOMET, STATE . SI ZE 

DIMENSION  ALPHA(2) ,X(SIZE) 

COMMON    ETA,K,PARAMG 

CALL    OVFLOW 
C 
C  PARAMETERS    FOR    NEGBIN 

K=DIST2K 

ETA=K/(DIST2M-1.0D0) 
C 
C    PARAMETER  FOR  GEOMET 

PARAMG=DLOG( 1 . ODO-1 .0  00/ DISTIM ) 
C 
C    STATE  ONE  TRANSITION  PROBABILITIES  FOR  MATRIX 

ALPHA(  1)=A1 

ALPHA(2)=1.0U0-A2 
C 
C    SELECT  INITIAL  STATE  FROM  STEADY-STATE  PROBABILITIES 

STATE=I 

C^LL  SRAND(IX,R, 1 ) 
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IF(R.LE.{(l.ODO-ALPHA(l))/(l.OD0-ALPHA(l)+/iLPHA(2)  ) 
C)  )  STATE=2 
C 
C    COMPUTE  'SIZE'  INTEREVENT  TIMES 

DO  2  1=1, SIZE 
C 
C    ENTER  MATRIX  AND  DETERMINE  TYPE  OF  NEXT  INTERVAL 

CALL  SRAND(IX,R,  1) 

IF(DBLE(R)  .LE.ALPHA(STATE) )  GO  TO  1 
C 
C    PICK  TYPE  TWO  VARIATE 

STATE=2 

X(I)=DFLOAT(NEGBIN(IX)) 

GC  TO  2 
C 
C    PICK  TYPE  ONE  VARIATE 

1  STATE=1 
X(I)=DFLOAT(GEOMET(IX)) 

2  CONTINUE 
RETURN 
END 


C  INTEGER  FUNCTION  GEOMET                   C 

C  C 

C    A.   PURPOSE:  C 

C  THIS  FUNCTION  GENERATES  VARIATES  FROM  THE           C 

C  GEOMETRIC  DISTRIBUTION  WHICH  IS  CF  THE  FCRM         C 

C  C 

C  X-1                         c 

C  M(X)=( 1-P)*{P)      ;0<P<l;X=l,2,...        C 

c  c 

C    B.   USAGE:  C 

C        THIS  FUNCTION  WAS  WRITTEN  TO  BE  USED  WITH  C 

C        SUBROUTINE  MODEL.  C 

c  c 

C        P=1-1/(1-DIST1M)  C 

C  C 

C        PARAMG=DLOG(P)  C 

c  c 

c    ^ , ^ , ...        c 

c  c 

INTEGER  FUNCTION  GE0MET*4  (IXJ 

IMPLICIT  REAL*8  (A-H,K,0-Z) 

REAL*4  R 

COMMON    ETA,K,PARAMG 
C 
C  CALCULATE    VARIATE 

CALL    SRAND(IX,R,1) 

RATIO  =  DLOG(DBLE(R))/PARAMG 

IRATIO  =  IDINT(RATIO} 
C 
C    ROUND  UP  IF  NON-INTEGER 

IF  (RATIC-DFLCAT( IRATIO). GT.O. ODD)  IRAT I  0= I  RAT  10+ 1 
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RETURN  VARIATE 
GEGMET=IRATIO 
RETURN 
END 
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C  INTEGER  FUNCTION  NEGBIN                   C 

C  C 

C  A.   PURPOSE:                                            C 

C  THIS  FUNCTION  GENERATES  VARIATES  FROM  THE  NEGATIVE  C 

C  BINOMIAL  DISTRIBUTION  WHICH  IS  OF  THE  FORM          C 

c  c 

C  (K-X-2)  K     X-1                         C 

C  M(X)=  (     )  *(1-P)  *(P)      ;K>0;0<P<l;x=l,2,...  C 

C  {  X'l  )  C 

C  C 

C    B.   USAGE:  C 

C  THIS  FUNCTION  WAS  WRITTEN  TO  BE  USED  WITH           C 

C  SUBROUTINE  MODEL.  C 

C  c 

C        P=1/(1+ETA)  C 

C  C 

c         , c 

c  c 

INTEGER  FUNCTION  NEGBIN*4  (IX) 

IMPLICIT  REAL*8  (A-H,K,G-Z} 

REAL*4  U(100),E,N0RM 

COMMON    ETA,K,PARAMG 

GAMMAD=O.ODO 

GAMMAI=O.ODO 
C 
C    DETERMINE  INTEGER  AND  DECIMAL  PARTS  OF  K 

IK=IDINT(K) 

DK=K-DFLOAT(IK) 
C 

C    CALCULATE,  IF  REQUIRED,  GAMMA  IK  VARIATE 
C    FROy  SUM  OF  IK  UNIT  EXPONENTIAL  VARIATES 

IFCK.LT.l.ODO)  GO  TO  9 

ET=0.000 

DO  8  M=l, IK 

CALL  SEXPON  (IX, E,l) 

ET=ET+DBLE(E) 

8  CONTINUE 
GAMMAI=ET/£TA 

C 

C    CALCULATE,  IF  REQUIRED,  GAMMA  DK  VARIATE 

C    USING  JOHNK'S  METHOD 

9  IF(DK.LE.O.ODO)  GO  TO  11 

10  CALL    SRAND    (IX, U, 2) 
UK1=DBLE(U(1) )*-( l.ODO/OK) 
UK2=DBLE(U(2)  J** ( 1 .000/ ( 1 .ODO-DK J  ) 
ZZ  =  UKH-UK2 

IF(ZZ.GE.l.ODO)     GO    TO    10 

CALL    SEXPON     (IX, E, 1) 

GAMMAD= (UK1/ZZ)*0BLE(EJ/ETA 
C 
C  TOTAL    GAMMA    VARIATE 

11  GAMMA=GAMMAO+GAMMAI 
IGAMMA=IDIM(GAMMA) 

IF( IGAMMA.GE.IOO)  GO  TO  50 
C 

C    CALCULATE  POISSGN  VARIATE,  DIRECTLY 
NN  =  0 

UT=1.0D0 

EGAMMA=DEXP(-GAMMA) 
20  CALL  SPAND  ( IX, U, 100) 
DO  30  M=l, 100 
UT=UT*DBLE(U(M) ) 
IF(UT.LE.EGAMMA)  GO  TO  40 
30  CONTINUE 
NN=NN+100 
GO  TO  20 
40  N  =  M-l  +  tJN 
GO  TO  60 
C 

C    CALCULATE  POISSON  VARIATE,  USING  NORMAL  APPROXIMATION 
50  CALL  SNCRM  (IX, NORM, 1) 
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N=IDINT(GAMMA-'0.2  5D0+(NGRM/2.0D0J**2+N0RM* 
C(0SQRT(GAMMA+0.12  5D0) )) 
C 

C    RETURN  VARIATE 
60  NEGBIN=N+1 
RETURN 
END 


C  SUBROUTINE  KSTEST  C 

c  c 

C  A.   PURPOSE:  C 

C        THIS  SUBROUTINE  CALCULATES  THE  PERIODOGRAM  OF  C 

C        INTERVAL  DATA,  ESTIMATES  THE  POWER  SPECTRAL  C 

C        DENSITY  PSD  AND  TESTS  THE  FIT  OF  THE  PERICDOGRAM  C 

C        TO  THE  PSD.  C 

C  c 

C  E.   USAGE:  C 

C        1-   ARGUMENTS:  C 

C  C 

C            X  -  INPUT  VECTOR  OF  INTERVAL  DATA  (REALMS)  C 

C  c 

C  IVEC  -  OUTPUT  VECTOR  OF  DESIRED  PERICDOGRAM     C 

C  POINTS  (REAL*8)  C 

C  c 

C  MEAN  -  OUTPUT  MEAN  OF  INTERVAL  DATA  (REAL*8)    C 

C  C 

C  VARIAN  -  OUTPUT  VARIANCE  OF  INTERVALS  (REAL*a)  C 

C  C 

C  XKSDN  -  OUTPUT  KOLMOGOROV-SMIRNQV  STATISTIC     C 

C  FOR  TEST  OF  FIT  (REAL*8)                 C 

C  C 

C  SIZE  -  INPUT  LENGTH  OF  VECTOR  X  (INTEGER)       C 

C  c 

C  MHAT  -  OUTPUT  ESTIMATE  OF  COVARIANCE  PARAMETER  C 

C  M  (REAL*8)  C 

C  C 

C  BHAT  -  PUTPUT  ESTIMATE  OF  COVARIANCE  PARAMETER  C 

C  BETA  (REAL*8)  C 

C  C 

C  2.   REQUIRED  SUBROUTINE:   FFTR  (IMSL  ROUTINE)  C 

C  C 

C  3.   caution:   VECTOR  X  OF  INTERVALS  IS  DESTROYED  C 

C  BY  FFTR.  C 

C  C 

c   , c 

c 

SUBROUTINE    KSTEST ( X , I VEC , ME  AN , V AR I  AN, KSDN , S  IZE , MHAT , 
CBHAT) 

ir'PLICIT    REAL*8     (A-H,0-Z) 

CCMPLEX*16  GAMN 

REAL*8  I  VEC,  ME  A.N,  MHAT,  I  B,  I  C  ,  KSDN  ,  KL,  KU 

INTEGER*4  SIZE, SS, HE 

DIMENSION  X(820) ,S(409) , I WK ( 249 5 ) , I VEC ( 409 ) 

DATA  PI/3. 14159265400/ 

GAMMA  1  =  0.0 DO 

GAMMA2=0.OD0 

KSDN^O.ODO 

MEAN=O.ODO 

VARIAN=O.ODO 

NN  =  SIZE  -  IDINT{DFL0AT(SIZE)/2)  -  1 

SS=  NN  -1 

HE  =  SIZE  -1 

JE  =  SIZE  -  2 
C 
C    CALCULATE  MEAN  AMD  VARIANCE  OF  INTERVALS 

DO  10  J=1,SIZE 

MEAN  =  MEAN  +  X( J  ) 
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VARIAN  =  VARIAN  +  X(J)**2 
10  CONTINUE 

MEAN  =  MEAN/SIZE 

VARIAN  =  (  VARIAN  -  SIZE  *  MEAN**2)  /  (SIZE  -  1) 
C 
C    CALCULATE  ESTIMATES  OF  M  AND  BETA 

DO  40  J=1,HE 

GAMMAl  =  GAMMAl  +  ( X( J  +  1) -MEAN  J  *  (X(J)-MEAN) 
40  CONTINUE 

DO  50  J=ltJE 

GAMMA2  =  GAMMA2  +  ( X( J+2 )-MEAN )  «  (X(J)-MEAN) 
50  CONTINUE 

GAMMAI=GAMMA1/SIZE 

GAMMA2=GAMMA2/SIZE 

BHAT  =  GAMMA2  /  GAMMAl 

MHAT  =  GAMMA1**2  /  GAMMA2 
C 
C    COMPUTE  FINITE  FOURIER  TRANSFORM  OF  INTERVAL  DATA 

CALL  FFTR(XtGAMN,SIZ£t IWK) 
C 
C    CALCULATE  PERIODOGRAM 

DC  20  J=3, SIZE, 2 

I=(J-l)/2 

IVEC(I)=(X(J)**2  +  X(J+1)**2)/(PI  *  SIZE) 
20  CONTINUE 
C 
C   TEST  PERIODOGRAM  FIT  TO  ESTIMATED  POKER  SPECTRAL  DENSITY 

OMEGA  =  DCQS(2  *  PI  /SIZE) 

PHAT=(VARIAK+2*MHAT*BHAT*(0MEGA-BHAT)/( l+BHAT**2-2* 
CBI-AT*OMEGA))/PI 

S(l)  =  IVEC{1)/PHAT 

DO  60  J=2,NN 

OMEGA  =  DC0S(2  *  PI  *  J  /  SIZE) 

PHAT=(VARIAK+2*MHAT*BHAT*(0MEGA-EhAT)/(l+BFAT**2-2* 
CBFAT*OMEGA))/PI 

S(J)=  S(J-l)  +  IVEC(J)/  PHAT 
60  CONTINUE 

DC  70  J=1,SS 

S(J)  =  S(J)/S(NN) 

KL=  DA6S(S(J)-(DFL0AT( J-1)/NN) ) 

KU=  DAdS(S(J)-(DFLDAT( J  )/NN) ) 

KSDN  =  DMAX1(KSDN,KL,KU) 
70  CONTINUE 

KSDN=KSDN*DSQRT(DFL0AT(NN) ) 
80  RETURN 

END 


C  SUBROUTINE  ESTIM8                       C 

c  c 

C         A         PURPOSE!  C 

C            *       THIS    SUBROUTINE    USES    THE    METHOD    OF    SCORING    TO  C 

C  STABILIZE    ESTIMATES    OF    THE    PARAMETERS     FOR    A                       C 

C  UNIVARIATE    TwO-STATE    SEMI-MARKOV    MODEL.                                  C 

c  c 

C    B.   usage:  C 

C  1.   ARGUMENTS:                                      C 

C  C 

C  M1,M2,S2,A1,A2  -  INPUT  INITIAL  ESTIMATES  FOR    C 

C  THE  MEAN  OF  TYPE  1  INTERVALS,  C 

C  MEAN  AND  STAr^DARO  CEVIATICN    C 

C  OF  TYPE  2  INTERVALS  AND  THE    C 

C  TRANSITICN  PROBABILITIES  -     C 

C  ALL  REAL«8                      C 

C  C 

C  SIZE  -  INPUT  NUMBER  OF  INTERVAL  DATA  POINTS     C 

C  (INTEGER)                                 C 

c  c 
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C  MEAN  -  INPUT  MEAN  OF  INTERV/iL  DATA  (REAL*8)  C 

C  C 

C  IVEC  -  INPUT  VECTOR  OF  PERIODOGRAM  POINTS  C 

C  (REAL*8)  C 

C  C 

C  ITERS  -  INPUT  NUMBER  OF  ITERATIONS  DESIRED  C 

C  PRIOR  TO  TERMINATION  IF  NO  CONVERGENCE  C 

C  (INTEGER)  C 

C  C 

C  CONVRG  -  CONVERGENCE  CRITERION  FOR  TERMINATION  C 

C  (REAL*8)   lO.OE-10  RECOMMENDED  C 

C  C 

C        2.   SUBROUTINES  REQUIRED:  C 

C  DMINV  (IBM  ROUTINE)  C 

C  DTERM  (IBM  ROUTINE)  C 

c  c 

C        3.   IF  NO  SCALING  IS  DESIRED,  SET  MEAN=1.0D0  .      C 

c  c 

C 

SUBROUTINE  ESTIM8 ( Ml ,M2, S2, Al, A2, SIZE, MEAN , I VEC, I TER8, 
CCCNVRG) 

I^'PLICIT  REAL*8  (A-Z) 

INTEGER*4  SIZE,NN,J,NE,JE,HE,I,L,M,N,ITE,ITER8 

DIMENSION  AVEC(5) ,AMAT(5,5) , MVEC ( 5 ) , MMAT (5,5) ,BVEC(5) 

DIMENSION  BMAT(5i5) ,L(5) ,LVEC(5) ,LMAT(5,5) ,CELT(5) 

DIMENSION  LLMAT(5,5),M(5) ,PVEC(5),IVEC(409) 

DATA  PI/3. 141592654D0/, N/5/ 

ITE=0 

NN=SIZE-I DINT (DFLOAT( SIZE) /2)-l 
C 

C    ZERC-OUT  VECTORS  AND  MATRICES 
90  DO  100  J=l,5 

AVEC(J)  =  0.000 

BVEC(J)  =  O.ODO 

LVEC(J)  =  O.ODO 

MVEC(J)  =  O.ODO 

DO  110  1=1,5 

AMAT( I, J)  =  O.ODO 

BMAT( I, J)  =  O.ODO 

LMAT(I,J)  =  O.ODO 

LLMAT( I,J)=  0.000 

MMATd  ,  J)  =  0.000 
110  CONTINUE 
100  CONTINUE 
C 
C    COUNT  ITERATIONS 

ITE  =  ITE  +  1 

WRITE  (6,345)  ITE 
C 
C    CALCULATE  ELEMENTS  OF  EQUATIONS 

CCl  =  1  -  Al 

CC2  =  1  -  A2 

CC3  =  CCl  +  CC2 

CC4  =  (Ml**2)  -  Ml 

CC5  =  Ml  -  M2 

CC6  =  CCl  *  CC2 

AE  =  ((CC2  *  CC4)  +  (CCl  *  S2**2))/CC3 

AVEC(1)=  MEAN  * ( ( 1-2*M 1 ) * (-CC2 ) ) /CC3 

AVEC(3)=  MEAN  *  (2*  32  *  CC1)/CC3 

AVEC(4J=  (CC2  *(CC4  -  52**2))/  (CC3**2) 

AVEC(5)=( (-CC1)*(CC4  -52**2))/  (CC3**2) 

AMAT(1,1)=  MEAN  *{2  *  CC2)/CC3 

AMAT(1,4)=  AVEC(l)  /  (CCS  *  MEAN) 

AMAT(1,5)=  (CCl  *(i-2*Ml) )/ (CC3**2) 

AMAT{3,3)=  MEAN  *(2*  CCl)/  CCS 

AMAT(3,4)=  (2  *  52  *  (-CC2))  /  CC3**2 

AMAT(3,5)=  AVEC(3)  /  (CC3  *  MEAN) 

AMAT(4,4)=  (2  *  AVEC(4))  /  CC3 

AMAT(4,5)=  ((A1-A2)  *  (CC4  -  52**2))/  CC3**3 

AMAT(5,5)=  (2  *  AVEC(5))/CC3 

NE  =  (CC6  «  CC5**2  )/CC3**2 
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MVEC(1)=  MEAN  *  (2  *  CC5  *  CC6)  /  CC3**2 

MVEC(2)=  -MVECd  ) 

MVEC(4) =(CC2*(A2-A1)*(CC5**2) )/CC3**3 

MVEC(5]=(CC1*(A1-A2)*(CC5**2) J /CC3**3 

MMAT(1,1)=MEAN*(2*CC6)/CC3**2 

MMAT(1,2)=-MMAT( 1, 1) 

MNAT(1,4)=(2*CC2*(A2-A1)*CC5)/CC3**3 

MMAT(2»4)=-MMAT( 1,4) 

MMAT(l,5)=(2*CCi*(Ai-A2)*CC5)/CC3**3 

MMAT(2,5)=-MMAT( 1,5) 

MMAT(2,2)=MEAN*(2*CC6)/CC3**2 

MMAT(4,4)=( ( (6*A2)+((-4)*A2**2)+(2*Al*A2)+( (-2)*AI)-2) 
C*(CC5**2) )/CC3**4 

MMAT(4,5)={( (2*( I-A1-A2+2*A1*A2) )-Al**2-A2**2 ) *CC5**2) 
C/CC3**4 

MMATC  5,5) =( {(6*A1)-4*A1**2  +  2*A1*A2-2*A2-2)*(CC5**2)) 
C/CC3**4 

DO  120  J=l,3 

DO  130  1=1,5 

AMAT( J,  I)=AMAT(J, I  )*MEAN 

MVAK J, I)=MMAT(J, I )*MEAN 
130  CONTINUE 
120  CONTINUE 

BETA=A1+A2-1.0D0 

LIKE=0.0D0 

DC  140  NE=1,NN 

COSINE  =  DCOS(  2  *  PI  *  NE/  SIZE) 

CC7  =  (1  +  BETA**2)  -  2  *  BETA  *  COSINE 

BE  =  1  +  (  2  *((BETA  *  COS  I NE) -BETA**2 )  ) /CC7 

BVEC(4)=  2*( ( (1+BETA**2)*CGSINE)-2*BETA)/CC7**2 

BVEC{5) =BVEC{4) 

BMAT(4,4)=2*(6*BETA**2-2*C0SINE*BETA**3+4*CCSINE**2-6* 
CBETA*C0SINE-2)/CC7**3 

BMAT(4,5)=BMAT{4,4) 

eMAT(5, 5)=BMAT(4,4) 
C 
C    CALCULATE  ESTIMATE  OF  POWER  SPECTRAL  DENSITY 

PE=  (AE  +  ME  *BE)/  PI 

CC8  =(IV£C(NE)  -PE  )/  PE**2 

CC9  =  (PE  -  2  *  IVEC(NE))  /  PE**3 
C 
C    CALCULATE  FIRST-ORDER  PARTIALS  OF  PSD 

DO  150  J=l,5 

PVEC(J)  =  (AVEC(J)  +  MVEC{J)*BE  +  BVEC(J)*  NE)  /  PI 
150  CONTINUE 
C 
C    CALCULATE  LIKELIHOOD  FUNCTION  VALUE 

LIKE  =  LIKE  -  (IVEC(NE)  /  PE  )  -  DLOG(PE) 
C 
C    CALCULATE  SCORES 

DO  160  JE=1,5 

LVEC{JE)=LVEC (JE)  +  CC8  *  PVEC(JE) 

DO  170  HE=JE,5 
C 
C    CALCULATE  SECOND-ORDER  PARTIALS  OF  PSD 

PJH=AMAT( JE,HE)  +  MMAT(JE,HE)  *  BE  +  MVEC(JE)  * 
CBVEC(HE)+MVEC(HE)*6VEC( JE)+BMAT( JE,HE)*ME 
C 
C    CALCULATE  SECGND-ORDER  PARTIALS  OF  LIKELIHOOD  FUNCTION 

LMAT( JE,HE)=LMAT(Jfc,HE )-CC9*PV EC ( J E ) *PVEC (HE ) -CC8* 
CPJh/PI 
C 
C    CALCULATE  EXPECTED  VALUE  OF  LMAT 

LLMAT( JE,HE)=LLMAT (JE,HE)+PVEC ( J E ) *PV EC ( HE ) /PE**2 
170  CONTINUE 
160  CONTINUE 
140  CONTINUE 
C 
C    FILL  IN  LOWER  TRIANGLE  OF  MATRICES 

DC  180  JE=1,4 

J=JE  +  1 

DC  190  HE=J,5 
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LMAT(HE,JE)=LMAT{ JE,HE) 
LLMAT(HE,JE)=LLMAT(JEtHE) 
190  CCNTINUE 
180  CONTINUE 

l^RITE(6j360) 
DO  200  1=1,5 

WRITE(6,370)  ( LM AT ( I , J ) , J=l , 5) 
00  185  J=l,5 
AMATd,  J)=-LMAT(  I, J) 
185  CCNTINUE 
200  CONTINUE 
C 

C    CALCULATE  INVERSE  AND  DETERMINANT  OF  LMAT 
CALL  DMINV(LMAT,N,D,L,M) 
WRITE(6,380)  D 
WRITE(6t390) 
00  210  1=1,5 

WRITE(6,370J  ( LMAT (I , J) , J=l , 5) 
210  CCNTINUE 
C 

C    CALCULATE  ADDITIVE  INCREMENTS  TO  ESTIMATES 
DO  240  1=1,5 
DELT(n=  O.ODO 
DC  250  J=l,5 

DELT(n=  DELT(I)  +  LMAT(I,J)  *  LVEC(J) 
250  CCNTINUE 
240  CCNTINUE 
C 
C    INCREMENT  ESTIMATES 

Ml=  Ml  +  OELT(l)  *  MEAN 
M2=  M2  +  DELT(2)  *  MEAN 
S2=  S2  +  DELT(3)  *  MEAN 
Al=  Al  +  DELT(4) 
A2=  A2  +  DELT(5) 
WRITE(6,430)  LIKE 
WRITE(6,440J  ( L VEC ( J )  ,  J  =  l , 5) 
WRITt(6,450)  M1,M2,S2,A1,A2 
C 
C    TEST  FCR  CONVERGENCE 

IF(DMAX1(DABS (DELT(l) ) , DABS ( DELT ( 2 ) ) , DABS ( CELT ( 3 ) ) , 
CDABS(DELT(4) ) ,DABS(DELT(5) ) ).LE.CCNVRG)  GC  TO  260 

IF(DMAX1(  DABS  (LVECd)  )  ,  DABS  (  LV  EC  (  2  )  >  ,  DABS  (  L  VEC  (  3  )  )  , 
CDABS(LV£C(4) ) ,DABS(LVEC{5) ) ) .LE.CONVRG)  GO  TO  260 
C 

C    TEST  FCR  NUMBER  OF  ITERATIONS 
255  IF(  ITE.LT. ITERS)  GO  TO  270 
280  WRITE  (6,455)  CONVRG 

GC  TO  265 
270  WRITE(6,460) 
C 

C    CONTINUE  IF  NECESSARY 
GO  TO  90 
260  WRITE(6,470) 
C 
C    TEST  LMAT  FOR  NEGATIVE  DEFINITENESS 

265  DC  275  1=1,5 
DO  267  JE=1,5 
DC  266  HE=1,5 

BMAT( Jc,HE)=AMAT( JE,HE) 

266  CCNTINUE 

267  CCNTINUE 

CALL  DTERMCI ,BMAT,D,N) 

AVEC( I ) =D 
275  CCNTINUE 

IFCAVEC (1) .LT.O.ODO. AND. AVEC(2) .GT.0.0D0.ANC.AVEC(3) 
C.LT.0.0D0.AND.AVEC(4i.GT.0.0D0.AND.AVEC(5) .LT. O.ODO) 
C  GC  TO  290 
285  WRITE(b,480) 

DC  295  1=1,5 

WRIT£(6,490J      I.AVECd) 
295    CONTINUE 

GC    TO    600 
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290  k^RITE(6,500) 
600  WRITE(6,400) 

00  220  1=1,5 

WBITE(6,370)  (LLMAT ( I , J ) , J= 1, 5 ) 

220  CONTINUE 
C 

C    INVERT  EXPECTED  VALUE  MATRIX 
CALL  DMINV(LLMAT,N,D,L,M) 
WRITE(6,410)  D 
WRITE(6t420) 
C 

C    TEST  INVERSE  FOR  POSITIVE  DEFINITENESS 
DC  230  I=lt5 

WRITE{6,370)  (LLMATd tJ) »J=lt5) 
DC  225  J=l,5 
AMAT( If J)=LLMAT(I , J) 
225  CONTINUE 
230  CONTINUE 

DO  610  1=1,5 
DO  620  JE=1,5 
DC  630  HE=1,5 
6MAT( JE,HE)=AMAT( JE,HE) 
630  CONTINUE 
620  CONTINUE 

CALL  DTERM(I,BMAT,D,N) 
AVECd)  =D 
610  CONTINUE 

IF(DMIN1(AVEC(1),AVEC(2),AVEC(3),AVEC(4) ,AVEC(5) ) .GT. 
CO.ODOJ  GO  TO  660 
640  WRITE( 6,510) 
DO  650  1=1,5 
WRITE(6,520)  I,AVEC(I) 
650  CONTINUE 

RETURN 
660  WRITE(6,530) 

RETURN 
345  FCRMATdHl  ,//////////,•  ITERATION  NUMBER', 14) 
360  FORMAT(//,'  NEGATIVE  MATRIX  OF  SECOND  PARTIALS',/) 
370  FCRMAT(/, 5020.10) 

380  FORMAT{//,'  DETERMINANT  OF  MATRIX ',//, D20  .10 ) 
390  FORMAT(//,«  INVERSE  MATRIX',/) 
400  FCRMAT(//,'  INFORM/^TICN  MATRIX',/) 
410  FORMAT!//,'  DETERMINANT  OF  INFORMATION  MATRIX',//, 

CD20.10) 
420  FORMATC//,'  INVERSE  INFORMATION  MATRIX',/) 
430  FORMATi//,'  LOG  LIKELIHOOD  FUNCTION  VALUE  =',D20.10) 
440  FORMATC//,'  SCORES  =',5D20.1G) 
450  FORMAT!//,'  PARAMETERS  =',5D20.10) 

455  FORMAT!//,'  NUMBER  OF  ITERATIONS  REQUESTED  WAS', 
C  REACHED  AND',//,'  CONVERGENCE  CRITERION  (  ', 
CD18.10,'  )  WAS  NOT  ACHIEVED.',//,'  RUN  TERMINATED') 
460  FORMAT!//,'  CONTINUING') 
470  FORMAT!//,'  COVERGENCE') 
480  FORMAT!  Ihi, /////, •  MATRIX  OF  SECOND  PARTIALS  IS  NOT  ', 

C'NEGATI VE  DEFINI TE' ,/) 
490  FORMAT!/,'  DETERMINANT  OF  PRINCIPAL  MINOR  SUBMATRIX  ', 
C'CF  ORDER', 12,'  OF  MATRIX  OF  SECOND  PARTIALS  =  ', 
CD20.10) 
500  FORMAT! IHl, /////, '  MATRIX  OF  SECOND  PARTIALS  IS  ', 

C'NEGATI VE  DEFINITE'  ) 
510  FORMAT!///,'  INVERSE  INFORMATION  MATRIX  IS  NOT  ', 

C  POSITIVE  DEFINITE' ,/ ) 
52C  FORMAT!/,'  DETERMINANT  OF  PRINCIPAL  MINOR  SLBMATRIX  ', 
C'OF  ORDER', 12,'  OF  INVERSE  INFORMATION  MATRIX  =  ', 
CD20.10) 
)  FORMAT!   .  . 
C  DEFINITE'  ) 
END 
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